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ABSTRACT 

We reconsider the pixel-based, "template" polarized foreground removal method within the context 
of a next-generation, low-noise, low-resolution (0.5 degree FWHM) space-borne experiment measuring 
the cosmological .B-mode polarization signal in the cosmic microwave background (CMB). This method 
was first applied to polarized data by the Wilkinson Microwave Anisotropy Probe ( WMAP) team and 
further studied by Efstathiou et al. We need at least 3 frequency channels: one is used for extracting 
the CMB signal, whereas the other two are used to estimate the spatial distribution of the polarized 
dust and synchrotron emission. No extra data from non-CMB experiments or models are used. We 
extract the tensor-to-scalar ratio (r) from simulated sky maps outside the standard polarization mask 
(P06) of WMAP consisting of CMB, noise (2 uK arcmin), and a foreground model, and find that, 
even for the simplest 3-frequency configuration with 60, 100, and 240 GHz, the residual bias in r is 
as small as Ar w 0.002. This bias is dominated by the residual synchrotron emission due to spatial 
variations of the synchrotron spectral index. With an extended mask with f s k y = 0.5, the bias is 
reduced further down to < 0.001. 

Subject headings: cosmic background radiation, cosmological parameters, early universe, inflation, 
gravitational waves 



1. INTRODUCTION 

Why study the .B-mode polarization of the cosmic mi- 
crowave background (CMB)? Detection of the primordial 
gravitational waves generated during inflation would give 
us a direct insight into the physical condition of the uni- 
verse when the energy scale w as close to the grand unifi- 
cation scale, ~ 10 16 GeV (see lLiddle fc Lvthll2009L for a 
recent review and references therein). While a direct 
detection of the primordial gravitational waves using, 
e.g., laser interferometers, seems not possible with the 
present-day technology, an indirec t detection using the 
5- mode polarization of the CM B (|Seliak fc Zaldarriaga! 
119971 iKamionkowski et"al1ll997t ) may be possible in the 
near future (most optimistically, within a few years), pro- 
vided that the energy scale of inflation at which the ob- 
served gravitational waves were generated was indeed as 
high as the grand unification scale. 

We often characterize the amplitude of gravitational 
waves (also known as tensor perturbations) using the so- 
called "tensor-to-scalar ratio," which is conventionally 
defined as 

_2(|/,+ P + |V x | 2 ) ! (1) 



(l^kl 2 ) ' 

where and h£ are the Fourier transform of the ampli- 
tudes of two linear polarization states of gravitational 
waves, and IZk is the primordial curvature perturba- 
tion, which is a scalar perturbation (hence the name, 
"tensor-to-scalar ratio"). It is TZ^ that seeded the ob- 
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served structure in the universe, as well as the dominant 
com ponent of the ob served CMB temperature anisotropy 
(see IWeinberl [20081 . for a recent review and references 
therein). 

The dominant, scalar part of the temperature 
anisotropy generates radial and tangenti al polarization 
patte rns around hot and cold spots ([Coulson et al.l 
119941 ). This is called the i?-mode polarization, and 
has been de t ected with high statistic al significance 
dBrown et al.l 120091: iChiang et all 120101 : lLarson et all 
120101 : IKomatsu et al]l20Tft IQUIET lf20inn . However, the 
-B-mode polarization, which cannot be generated by the 
scalar perturbations but can be generated by the ten- 
sor perturbations, has not been found yet. The current 
95% upper limit on the tensor-to-scalar ratio is r < 0.24, 
which mainly comes from the upper limit on the ten- 
sor contributio n to the temperature anisotropy on large 
angular scales (jKomatsu et al.ll2010f ). 

Given the upper limit on r, one can calculate the ex- 
pected level of the .B-mode power spectrum (see Fig- 
ure [J). For r = 0.24, the -B-mode power spectrum is 
smaller than the .E-mode power spectrum by a factor of 
10 at the first bump (created by electrons at z < 10). 
At the second bump (created by electrons at z ~ 1090), 
the .B-mode power spectrum is smaller than the -E-mode 
power spectrum by a factor of 50. It is the smallness of 
the .B-mode power spectrum that makes the detection of 
this signal challenging. 

There are three sources of noise for .B-mode detec- 
tion: (1) Detector noise; (2) Galactic foreground emis- 
sion; and (3) Gravitational lensing. In this paper, we 
shall focus on the Galactic foreground. We use a map- 
based method for reducing the Galactic foreground, and 
study how the residual foreground limits a measurement 
of the primordial B-mode polarization. The foreground 
reduction technique we use is motivated by the "template 
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Fig. 1. — _E-mode and B-mode polarization power spectra. The 
diamonds, triangles, stars, a nd squares show the WMAP seven- 
year data IILarson et al. 201j]), the QUaD final data l|Brown et alj 
l2009t) . the BICEP two -year data llChiang et al.l [20101) , and the 
QUIET 43 GHz data iQU IET 20105, respectively. The upper 
solid line shows the scalar B-mode power spectrum of the WMAP 
seven-year best-fit model. The dashed lines show the primordial 
B-mode power spectra with the tensor-to-scalar r atio of r = 0.24, 
which corresponds to the current 95% upper limit IIKomatsu et alj 
[20101 ). as well as of r = 0.03 and 0.003. These lines are linearly pro- 
portional to r. The dotted line shows the secondary B-mode power 
spectrum expected to be generated by the w eak gravitational lens- 
ing e ffect converting E modes to B modes (Zaldarriaga <fc Scliak 
119981) . This line is fixed (by the WMAP seven-year best-fit model) 
and acts as noise for the primordial B-mode detection. The lens- 
ing contribution becomes comparable to the primordial bump at 
I = 10 and 100 for r = 0.003 and 0.03, respectively. 



clean i ng method" used by th e WMAP team (|Page et al.1 
[20071 IGold et al.|[2009l l2010lh T his method was further 
investigated bv lEfstathiou et alj ([2OO90 in the context of 
the Planck mission. We shall study this technique in the 
context of a next-generation, low-noise, low-resolution 
(0.5 degree FWHM) space-borne experiment. 

There is a large body of literature on the issue of polar- 
ized foreground cleaning for the i?-mode detection. Our 
method is one specific (and relatively simpler) example. 
For t he other methods in the literature, see r eview arti- 
cles (jDunklev et al.l 120091 ; iFraisse et al.ll2008ft and refer- 
ences therein. 

This paper is organized as follows. In Section [21 we 
show how the detector noise and the lensing noise influ- 
ence the statistical errors on r. In Section [31 we describe 
our method for estimating r in the presence of the Galac- 
tic foreground and the dominant scalar i?-mode polariza- 
tion. In Section 21 we describe our simulation including 
CMB, detector noise, and foreground. In Section [51 we 
present the main results of this paper. We conclude in 
Section [5] 

2. DETECTOR NOISE AND LENSING NOISE 

Before we study the effect of the foreground, we show 
how the detector noise and the lensing noise influence 
our ability to detect r. The detector noise enters into 
the likelihood of r via the noise power spectrum, N^ B . 
Assuming white noise, we write the noise power spectrum 
as 



N BB 
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10800 /iK arcmin 



/zK str, 



(2) 



Fig. 2. — B-mode polarization signal and noise power spectra. 
The dashed lines show the primordial B-mode power spectra with 
the tensor-to-scalar ratio of r = 0.01 and 0.001, while the dotted 
line shows the secondary B-mode power spectrum from gravita- 
tional lensing. We also show the noise power spectra (Equation l[2]|) 
— 1/2 

for w p = 2 and 10 /iK arcmin with a Gaussian beam window 
function of 6> FW HM = 30 arcmin, i.e., JV^e'^FWHM/t 81 " 12 ). 



— 1/2 

where w p is the noise in Stokes parameters Q or U per 
pixel whose solid angle, f2 p i x , gives \/£l p i x = 1 arcmin. 
This quantity is useful because one can compare various 
experiments on the same scale. 

Current and future experiments use many (of order 
10 3 — 10 4 ) detectors to reduce the noise equivalent tem- 
perature (NET) down to a few uK arcmin level. Is 
this sufficient for detecting primordial B modes? For 
comparison, the expected sensitivity of Planck com- 
bining 70, 100, and 143 GHz is w~ 1/2 = 63 /iK ar- 
cmin (see, e.g., Append ix A of iZaldarriaga et all 120081 : 
IPlanck Blue Book 1 127)051 1 . 

In Figure [21 we compare the noise power spectra for 

— 1/2 

w p =2 and 10 /iK arcmin to the primordial and lens- 
ing B modes. For r = 10 -3 and the 10 /iK arcmin noise, 
only a few modes {I = 2, 3, and 4) are above noise. For 
the 2 /iK arcmin noise, the noise power spectrum is be- 
low the lensing i?-mode power spectrum, and thus noise 
is no longer the limiting factor (unless we "de-lens" maps 
and remove the lensing noise) . How would this influence 
our ability to detect r? 

To see this, let us calculate the likelihood of r for a 
given noise level. For simplicity, we assume that we 
cover the full sky and the noise per pixel is homoge- 
neous Q Then, one can write down the probability dis- 
tribution function of the measured i?-mode power spec- 
tr um, CP B , for a given value of r as (e.g., Equation (8) 
of lHamimeche fc L ewis 20Ql) 



-2 In P(Cf B \r) 



:(2Z + 1) 
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ln(rcf ^ + of 
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cf + 



N BB 



1 



21 + 1 



HC? B ) 



(3) 



3 We assume this only for producing Figures [3] and [4] For the 
main analysis, we include inhomogeneous noise, foreground, and a 
partial sky coverage. 



Simple foreground cleaning algorithm for detecting i?-mode polarization of CMB 



3 



.^fV, mox =2 

/ / with lensing 

; h w/o lensing 


Lox = 5 

; r ft '-. lOjtiKorcmin 
; A \ 
; A \ '**. 
-" A 2/iKarcmin\ \ 




' mo " =100 " 



1.0 
0.8 
I 0.6 

> 0.4 

0.2 

0.0 
1.0 

0.8 

I 0.6 

> 0.4 
0.2 
0.0 

0.0 0.5 1.0 1.5 2.0 2.5 0.0 0.5 1.0 1.5 2.0 2.5 
Tensor-to-scalar Ratio, r [xlO 3 ] 

Fig. 3. — Effect of detector noise and gravitational lensing on the 
likelihood for r (foreground is not included). The input value of 
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multiplied by 10 3 . In the top left, top right, bottom left, and bot- 
tom right panels, we sum the multipoles up to i max = 2, 5, 10, and 
100, respectively. For the detector noise level, we use 2 fiK arcmin 
for the solid and dashed lines, and 10 fiK arcmin for the dotted 
lines. The solid and dotted lines include the gravitational lensing 
contribution to the total noise, while the dashed lines do not. Even 
if we set the detector noise to be zero, the solid lines do not change 
very much: for r = 10~ 3 , the gravitational lensing effect prevents 
us from measuring r accurately beyond I ~ 10. Note that a single 
multipolc, I = 2, is enough for us to detect r = 10~ 3 if the detector 
noise is smaller than 10 fiK arcmin. 
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Fig. 4. — Fractional error, cr r /r; nput , on the determination of 
the value of r as a function of maximum multipoles, i max . Here, 
ay is the square-root of the second-order moment of the likelihood 
function given by Equation JB). (The foreground is not included. 
The full sky coverage is assumed.) From the top to bottom lines, 
we show ri n p U t = 0.001 with and without the lensing noise, and 
Input = 0.01 with and without the lensing noise, respectively. For 
the instrumental noise level, we have used 2 fi¥L arcmin. Note 



that oy/nnput ~ 1 at l„ 



3 does not mean that we do not 



where cf w is the primordial i?-mode power spectrum 
from gravitational waves with r — 1, and cf is the sec- 
ondary B mode from gravitational lensing. We then 
use Bayes' theorem to calculate the likelihood for r as 
C(r\CP B ) oc P{C BB \r). To calculate the likelihood, 
we set the measured power spectrum to 

^inputC^ ~h ( 

multipoles up to l„ 



detect r; on the contrary, we detect r with high significance even 
at Z m ax = 2. Rather, it just means that the likelihood for r is highly 
non-Gaussian and has a long tail toward large values of r (see the 
top left panel of Figure [3). In other words, we detect r with high 
statistical significance, but the value of r is not determined very 
well. 



Nf B t 



be C BB = 
and sum the log-likelihood over 



lihood: 



drC(r)r 2 — 



drC(r) 



(5) 



ln£(r) = Vln£(r|Cf s ). 



1=2 



(4) 



Here, we have assumed that the likelihood is normalized 
such that J Q drC{r) = 1. One should be careful when 
interpreting this quantity. For / max = 2, ay would be 



Figure [3] shows the likelihood of r for the input value greater than the input value, 



of r 



input 



10" 3 and L 



2, 5, 10, and 100. One useful 



input 



= 10 3 ; however, this 



number to keep in mind is that a single multipole, I = 2, 
is sufficient for detecting r = 10 ~ 3 , if the noise is smaller 
than 10 fiK arcmin. However, the precision on r does not 
improve beyond 1 = 5. This is apparent also in Figure (2) 
the noise power spectrum exceeds the signal at I > 5. 

We can improve the precision further if we lower the 
noise level to, say, 2 fiK arcmin. Even so, the gravita- 
tional lensing prevents us from improving on the preci- 
sion beyond I <~ 10 if r = 10~ 3 . (If there were no lensing 
in the universe, we would be able to continue to improve 
on the precision, as indicated by the dashed lines.) In 
fact, 2 fiK arcmin is essentially the same as zero detector 
noise, as the lensing term dominates the error budget. 
Again, this is apparent in Figure [21 

Of course, these results are overly optimistic, as the er- 
ror would be dominated by the foreground rather than by 
the detector noise. Nevertheless, it is still useful to know 
what would be possible when we ignore the foreground. 

To quantify the precision on r, it is convenient to use 
the variance, of , given by the second moment of the like- 



does not mean that we cannot detect r. This just means 
that the distribution is highly non-Gaussian and has a 
long tail toward large values of r (see the top left panel 
of Figure [3]). For large values of l maXi e.g., l max > 10, 
the distribution of r becomes approximately a Gaussian, 
and thus the value of ay may be interpreted as the size 
of the usual la error bar. 

Figure [4] shows the fractional error, a r /ri nput , on the 
determination of the value of r as a function of i max - 
First, as one may expect from Figure El the fractional 
error for ri npu t = 10 -3 saturates at l max ~ 10 and docs 
not improve further due to the lensing noise. For this 
case, while we can detect r with high statistical signif- 
icance, we can determine the actual value of r to only 
^ 20%. For rinput = 10 -2 , we can determine the value 
of r to ~ 4% at Z max = 200 (beyond which the fractional 
error no longer improves due to the lensing noise). 

This study gives us an estimate of statistical errors 
on the measured values of r. On the other hand, the 
Galactic foreground gives us systematic errors (and bias). 
Now we shall turn to the foreground issue, which is the 
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main subject of this paper. 

3. PIXEL-BASED FOREGROUND REMOVAL METHOD 
3.1. Motivation 

The basic idea behind our methodology is simple: we 
have (at least) 3 polarized components on the sky that 
we know and have been detected by the WMAP: CMB, 
synchrotron emission, and thermal dust emission. As the 
synchrotron dominates at lower frequencies and the dust 
at higher frequencies, we use one map at a low frequency 
and another map at a high frequency as the foreground 
"templates." We put the quotation marks here because 
these maps also contain the CMB. No external template 
maps are used in our method. 

The WMAP team has applied this method for mod- 
eling the synchrotron: they used the lowest frequency 
(K-band, 23 GHz) map as a template, fitted it to the 
higher frequency maps (Ka, Q, V, and W bands), and 
subtracted from those maps. One can write this opera- 
tion as 

[Q , j v n, v) = [Q, £/](!/) -asMQ, = 23 GHz) 

l-as{y) 

where Q' and U' are the template-cleaned Stokes Q and 
U maps, respectively, and as is the best-fit synchrotron 
coefficient for a given frequency v. The denominator ac- 
counts for the fact that the K-band map also contains 
the CMB signal. 

However, the WMAP team had to rely on an external 
map for modeling the dust emission, as the highest fre- 
quency, the W band (94 GHz), was not high enough for 
being a good template of the polarized dust emission. 
This issue would probably be resolved by the Planck 
satellite, which ha s higher frequen c y cha nnels such as 
217 and 353 GHz. lEfstathiou et al.l (|2009t) have studied 
this by using a simulated Planck 217 GHz or 353 GHz 
map as a template for dust, and a simulated 30 GHz 
map as a template for synchrotron. They find that this 
simple method removes the foreground efficiently, bring- 
ing the bias in r down to a few times 10 -3 , which is 
much smaller than the expected statistical uncertainty 
on r from Planck, oy = 0(1O~ 2 ). 

The goal of this paper is to put this method in the 
context of a next-generation, low noise (2 fiK arcmin) 
polarization satellite experiment, and see if this method 
yields a promising result for measuring r ~ 10~ 3 (which 
is easy to detect in the absence of foreground, as we just 
saw in Section [2J . 

3.2. "Template" cleaning method 

Our methodology is sim ilar to that given in Section 4.2 
of lEfstathiou et al'j (pOOl . 

The main parameter that we wish to extract from data 
is the tensor-to-scalar ratio, r. (We do not vary the ten- 
sor tilt, n t .) The foreground coefficients, a, are nuisance 
parameters that we wish to marginalize over. The fore- 
ground coefficients may be spatially varying. 

Another nuisance parameter (for detecting B modes) 
is the amplitude of the scalar S-mode power spectrum, 
which is by far the dominant source of CMB polarization. 
The signal power spectra are thus given as 



tensor. BB 



(8) 



where q denotes the power spectra with s = 1 and r = 1. 
The fiducial value of s is s = 1. 

We shall maximize the following likelihood function for 
estimating r, s, and o^: 



C(r,s,ai) oc 
where 



exp 



\x'(ai) T C 1 (r,s,a i )x'(a i ) 



V\c{ 



r.s.cti 



luKv)-Ei^)lQ,u}(vT mplate ) 



(9) 



(10) 



is a template-cleaned map. This is a generalization of 
Equation ([6]) for a multi-component case. In this paper, 
i takes on "S" and "D" for synchrotron and dust, respec- 
tively, unless noted otherwise. For definiteness, we shall 
choose: 



V = 
template 



template 
D 



100 GHz, 
= 60 GHz, 
: 240 GHz. 



EE 



■ SC. 
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,EE , tensor, EE 
+ TC, , 



(7) 



These choices are somewhat arbitrary, but our prelim- 
inary optimization study indicates that this is a good 
configuration for achieving a smaller bias in r. A fuller 
optimization study, including more frequency channels, 
would require a more detailed specification of a given ex- 
periment (e.g., how many detectors one can fit in a given 
focal place; how low the detector noise can be as a func- 
tion of frequencies), which is beyond the scope of this 
paper, but will be presented elsewhere. 

The covariance matrix in pixel space, C, for Stokes Q 
and U maps is given as 

C(r, s, ai ) = rc tensor + sc scaIar + ^1 + ^2 (n) 

(1 - >J 4 ai) 2 

where c is the signal covariance matrix calculated from 
the theoretical power spectra, q, (see Appendix |A"]) and 
the noise matrices, N\ and N2, are a noise covariance of 
a smoothed map (which is not diagonal) before the tem- 
plate cleaning is applied, and a small artificial diagonal 
noise matrix for a matrix regularization, respectively (sec 
Section H3] for details). 

For simplicity and clarity, wc have ignored noise in 
template maps. For, if wc assume that all three channels 
are similar in detector noise level, it is a good approxi- 
mation, as «d ~ 0.08 and as ~ 0.25, and the fractional 
contribution of the template noise to the covariance ma- 
trix is given by a 2 , i.e., 6% effect in the derived error 
bars. Not e that this is equivalent to ignoring P in Sec- 
tion 4.2 of lEfstathiou et al.l (|2009]) . 

4. SIMULATION 

4.1. CMB and detector noise 

For CMB, we first generate the scalar and t ensor polar- 
izatio n power spectra using the CAMB code dLewis et al.l 
l2000f ) with and without lensing contributions. We then 
generate Stokes Q and U maps at the Healpix resolu- 
tion of iVgide = 128. The signal map has been smoothed 
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Fig. 5. — Foreground maps from the Planck Sky Model (PSM; vl.6.2). The top left and top right panels show the polarization intensity 
maps (P = Q 2 + U 2 in units of fiK) of synchrotron and dust, respectively. The dust polarization intensity has been multiplied by a 
factor of three to better approximate a more recent version of PSM. The lower left and lower right panels show the synchrotron index /3s 
and the dust index , respectively. Note a small range shown for j3 n : the dust index does not vary much, but this is a built-in assumption 
of the PSM vl.6.2. 



with a 30' beam (FWHM), representing a low- angular- 
resolution CMB polarization satellite experiment target- 
ing the primordial B modes. 

To this smoothed signal map, we add random Gaussian 
noise given by <Jq / \J N t, s (n) per pixel in the direction of 

— 1/2 

n. Here, ctq is related to noise w p as 



-1/2 



00 



10800 fiK arcmin Jn^N^E. K^m) 



(12) 

where 7V pix = 12(128) 2 = 196608 is the total number of 
pixels at A^dc = 128, and -ZV b s is the number of observa- 
tions per pixel. We adopt -/V n b s from t he "EPIC low-cost" 
(EPIC-LC) design (|Bock et al.ll2008[ ). The noise is high- 
est on the ecliptic plane and lowest on the ecliptic poles, 
similar to the iV bs pattern of the WMAP. Note that the 
absolute value of Aob s will cancel out in er / V N h s (h) if 
we use the above formula: only the spatial distribution is 
taken from A^bs, and the overall noise level is set by the 



— 1/2 

assumed value of w p . We shall use w p 1 '" =2 fiK ar- 
cmin for the rest of this paper. For this low noise config- 
uration, the results are not sensitive to the details of the 
N hs pattern. 

As we described at the end of Section [3l noise in tem- 
plate maps (at = 60 GHz and v& = 240 GHz) makes 



-1/2 



only a small contribution to the final covariance matrix. 
Therefore, for simplicity we add noise only to our CMB 
channel at 100 GHzQ 

We then apply an additional Gaussian smoothing to 
this signal-plus- noise map with 9.16 degrees (FWHM), 
which is 2.5 times the pixel size at A^idc = 16, and re- 
sample the smoothed map to A^ido = 16. Finally, as 
the smoothed map at A^idc = 16 is dominated by the 
scalar i5-mode signal at all angular scales supported by 
the map resolution, the covariance matrix of this map 
is singular. In order to regularize the covariance ma- 
trix, we add an artificial, homogeneous white noise of 
0.2 fiK arcmin such that the map becomes noise domi- 
nated at the Nyquist frequency, i max = SA^do — 1 = 47. 

4.2. Foreground: Planck Sky Model 

For the Galactic foreground model, we use the Planck 
Sky Model (PSM; vl.6.2) developed by the Planck Com- 
ponent Separation Working Group (Working Group 2). 
iLeach et al.l (120081) descr ibe the PSM for temperature, 
and lDunklev et aLl (|2009t) for polarization. 

The polarized synchrotron and dust emission are mod- 

4 Note that noise in templates cannot be ignored when we try 
to find an optimal combination of 3 frequencies. We ignore noise 
in templates here because we have done our preliminary optimiza- 
tion already. A fuller exploration of template noise along with the 
frequency optimization will be given elsewhere. 
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Fig. 6. — B-modc polarization signal and foreground power spec- 
tra. The dashed lines show the primordial i?-mode power spec- 
tra with the tensor-to-scalar ratio of r = 0.01 and 0.001, while 
the dotted line shows the secondary B-mode power spectrum from 
gravitational lensing. We also show the B-mode power spectra of 
the dust (dot-dashed line) and synchrotron (triple-dot-dashcd line) 
emission at 100 GHz outside the WMAP P06 mask (N sidB = 128). 
(The measured foreground power spectra have been divided by 
/sky = 0.733 to approximately correct for the mask, as well as by 
the pixel window function at N s i de = 128.) The total (synch+dust) 
power spectrum is measured from the total map, which is slightly 
larger than the sum of the synchrotron and dust power spectra, 
as these two foreground components are spatially correlated. Note 
that the original PSM vl.6.2 dust map has the average polarization 
fraction of 1.5%, but we have multiplied the dust map by a factor 
of 3 to approximate a more recent dust template map adopted by 
the Planck collaboration. 

eled as power-laws in antenna temperature: 
[Q sy „ch, U syach ](u, fi)=g(v)[Q^ h , U^ h ](S0 GHz, h), 

x (sogs) ' (13) 

[Qdust, U iwt ](v, n)=g(v)[Q2™, ^1(94 GHz, n), 

(14) 



94 GHz 



Here, Q PSM and 77 PSM arc the PSM Stokes parame- 
ters in units of antenna temperature, and g(i>) = (e x — 
l) 2 /(x 2 e x ) where x = hv/k B T CM B = ^/56.780 GHz con- 
verts the antenna temperature to thermodynamic tem- 
perature. (Q and U are in units of thermodynamic tem- 
perature.) 

For synchrotron, the position-dependent spectral in- 
dex, (3s (n) , is calcul ated from the Haslam 408 MHz map 
(|Haslam et al.lll98l"l ) a nd the three-year WMAP temper- 
ature map at 23 GHz ([Page et al.ll2007|). The template 
maps at 30 GHz are taken from iMiville-D cschcnc s et al.l 
(l2008l) . 

For dust, the position-dependent spectral index, 
/3d (w) , as well as t he unpolarized intensity map are taken 
from Model 8 of iFinkbeiner et al.l (|1999|) . The polar- 
ization angles of dust approximately follow those of the 
synchrotron maps. The original PSM dust map has the 
average polarization fraction of 1.5% over the full sky, 
but we will multiply this map by a factor of 3 to ap- 
proximate a more recent dust map used by the Planck 
collaboration. 

Top panels of Figure [5] show the amplitude (P = 
\/Q 2 + U 2 ) of polarization intensity of synchrotron and 
dust at 100 GHz, while the bottom panels show the 



spectral indices, /3g and /3d- After adding the above 
foreground maps (smoothed with a 9.16-degree beam at 
-^Vsidc = 128 and degraded to iV s i d o=i6) to the CMB-plus- 
noise map, we mask th e simulated sky by the WMAP 
P06 mask (f sky = 73%) jPaje et al.ll27)07h . 

The norm of the pixel vector, [Q,E/], is 2259 x 2, where 
2259 is the number of pixels outside the P06 mask. In 
order to mask the covariance m atrix, we us e the tech- 
nique described in Appendix D of lPage et al.l (|2007|) : we 
compute an inverse of 6144 x 6144 matrix a nd reduce it 
to 451 8 x 4518 matrix using Equation (D7) of lPage et al.l 
(|2007D . (Note that there is a typo in this equation: D 
should be replaced by D~ l .) 

In Figure [6l we show the .B-mode power spectra mea- 
sured from the PSM (N sidc = 128) at 100 GHz outside 
the P06 mask. The total foreground power spectrum has 
1(1 + l)C l BB /{2n) w 10" 1 fiK 2 at I < 10, which is 250 
and 2500 times larger than the primordial B-mode spec- 
tra with r = 0.01 and 0.001, respectively. The problem 
seems formidable; however, as we show below, the sim- 
ple cleaning method can reduce the foreground-induced 
bias in rto Arw 0.002(< 0.001) with the P06(extended) 
mask. 



5. RESULTS 



5.1. 



Fixing the scalar E-mode amplitude 

Before we use our full likelihood function given by 
Equation ([9]), let us first try a simpler version and show 
that it actually fails. 

For the moment (only within this subsection), we fix 
the amplitude of the scalar E modes, i.e., s = 1, and 
consider cleaning dust using a map at 240 GHz. (Syn- 
chrotron will not be discussed in this subsection.) Our 
model is thus 



[Q,U] (100) = CMB 
[Q, U] (240) = CMB ■ 



Dust(100) 
Dust (240). 



Noise, (15) 
(16) 



As we described at the end of Section [3l we ignore noise 
at 240 GHz. We then fit the 240 GHz map to the 100 GHz 
map: 

[<?', 1/1(100) = [Q, C/](100) - ao[Q, 17] (240). (17) 

Minimizing x 2 = [Q', U'] T C- 1 [Q', U'] with respect to 
«d gives the following least-square solution: 



[Q,t/] T (100)C- 1 [Q,t/](240) 
[Q,[/]^(240)C- 1 [Q,C/](240)- 



(18) 



As the polarization signal is dominated by scalar E 
modes, we can set r = when computing the covariance 
matrix C in this equation. (In practice, we used Hnput-) 
Finally, we maximize the likelihood given in Equation ([9]) 
with respect to r, with s = 1 and qd given by the above 
least-square solution. 

The left panel of Figure [7] shows the values of r and 
ao obtained from many random realizations of noise 
and CMB skies. (The input tensor-to-scalar ratio is 
Hnput = 0.003.) There is a clear correlation between r 
and ao , indicating a failure of this algorithm. This cor- 
relation is caused by a chance correlation between fore- 
ground and the dominant sc alar E modes (|Chiang et al.l 
2008; Efstathiou et al. 2009). The correlation disappears 
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Fig. 7. — Correlation between the tensor-to-scalar ratio, r, and the amplitude of dust, ojj. Left: The amplitude of the scalar E modes 
is held fixed at s = 1. Right: The amplitude of the scalar E modes is treated as a nuisance parameter and marginalized over. The input 
tensor-to-scalar ratio is ri nput = 0.003 and, for this figure only, the original PSM dust map (with an average polarization of ~ 1.5%) is 
used. 



TABLE 1 

Dust-only Test 



^input 


mean(r) b 


std(r) c 


0.001 


0.0011 


0.0003 


0.003 


0.0030 


0.0005 


0.010 


0.0102 


0.0010 


0.030 


0.0296 


0.0021 


0.100 


0.0991 


0.0057 



a Input values of the scalar-to- 
tensor ratio for simulations (64 
realizations for each r; nput ). 
b Mean of the recovered maxi- 
mum likelihood values of r. 
c Standard deviation of the 
recovered maximum likelihood 
values of r. 



when we set Cf E = 0. This result motivates our treating 
the amplitude of scalar modes as a nuisance parameter. 

The right panel of Figure [7] shows the results when 
s is treated as a nuisance parameter and marginalized 
over. For this, we have maximized the likelihood given 
by Equation (|9]) by varying r, s, and cud simultaneously. 
The correlation between r and qd has disappeared. 

How well was dust cleaned? We have repeated this one- 
component foreground cleaning test for various values of 
Hnput from 0.001 to 0.1. The results are shown in Table [TJ 
in all cases, the method recovers r successfully. 

5.2. Cleaning synchrotron in multi-region 
We are now ready to include synchrotron. Our model 

is 

[Q, U] (60) = CMB + Synch(60) + Dust(60) (19) 
[Q, U] (100) = CMB + Synch(lOO) + Dust(100) 

+Noise (20) 
[Q, U] (240) = CMB + Synch(240) + Dust(240) (21) 



It turns out cleaning synchrotron is more challenging 
than cleaning dust, as the spatial distribution of syn- 
chrotron tends to be more extended above the Galactic 
plane than that of dust (see the top panels of Figure [5]) . 
We start by adding a mock synchrotron model (MSM) 
map to the PSM dust map. The MSM map has the same 
synchrotron polarization intensity across the sky as PSM 
at 30 GHz, but has a spatially invariant spectral index 
of P = —3.0 (where (3 = —3.0 is the average of spatially 
varying spectral index of PSM). With MSM plus PSM 
dust, r is recoverd successfully; mean(r) = 0.0012 and 
0.0031 for ri npu t — 0.001 and 0.003. (see the second and 
third columns of Table [2]) 

Even more problematic is the spatial variation of the 
synchrotron spectral index (see the bottom left panel of 
Figure [5]) , which causes a mismatch between a template 
map at 60 GHz and the actual synchrotron distribution 
at 100 GHz. When we use a single synchrotron coef- 
ficient, as, for the whole sky for the PSM model even 
without dust (= synchrotron only), we find a bias in r 
of order Ar w 0.002: mean(r) = 0.0028 and 0.0120 for 
Input = 0.001 and 0.01, respectively (see the fourth and 
fifth columns of Table E]). 

One way to mitigate this issue woul d be to extend the 
Galactic mask (Efstath iou et al.ll2009D . In addition, one 
may give up using a single synchrotron amplitude for the 
whole sky, and use multiple amplitudes depending on the 
locations on the sky0 In this paper, we 

(Method I) continue to use the P06 mask, but di- 
vide the sky using the Healpix map with 
TVgide = 2, as shown in Figure and 

(Method II) extend the mask to f s k y = 50% and di- 
vide the sky using the Healpix map with 

5 Ultimately, the best way to mitigate this issue would be to 
obtain and use information on the spatial distribution of the syn- 
chrotron spectral index. 
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TABLE 2 

MSM and Synchrotron-only Tests (global and 48 regions) 



'"input 


MSM 


b 


Global 


48 Regions 


mean(r) d 


std(r) c 


mean(r) d std(r) c 


mean(r) d std(r)° 


0.001 


0.0012 


0.0004 


0.0028 0.0005 


0.0024 0.0005 


0.003 


0.0031 


0.0006 


0.0049 0.0008 


0.0046 0.0007 


0.010 






0.0120 0.0011 


0.0115 0.0011 



a Input values of the scalar-to-tensor ratio for simulations (64 realizations for 

each r input ). 

b MSM plus PSM dust. 

c PSM synchrotron only. 

d Mean of the recovered maximum likelihood values of r. 

c Standard deviation of the recovered maximum likelihood values of r. 



48 a s regions with the P06 mask (/ aiJ -Wo] far Method I 




I 



12 n s regions with extended mask (/^-50%) for Method II 




Fig. 8.— (a) left: Method I. The P06 masked sky has been 
divided into 48 regions based on the Healpix map with N s id c = 2. 
fsky of the mask is 73%. (b) right: Method II. The extended- 
masked sky has been divided into 12 regions based on the Healpix 
map with N a i^ e = 1. f a k y of the mask is 50%. 



iVs 



de 



1, as shown in Figure [SJa. 



We give the details of our definition of the extended mask 
in Appendix [B] In short, we choose the threshold polar- 
ization intensity values at 60 and 240 GHz above which 
the pixels arc masked, such that we retain 50% of the 
sky. 

While we would probably do a better job at cleaning 
synchrotron if we divide the sky according to our knowl- 
edge of the polarized synchrotron measured by WMAP, 
for this paper we prefer to explore a simpler algorithm 



and see how well we can recover r. 
Each region I will be cleaned as (c.f., Equation (flUl) ) 



x, = 



[Qj, U>](100) - ap[Qi, t/ 7 ](240) - o|[Qj, lfr](60) 

1 — q;d — Oig 

(22) 

Note that we still use a single amplitude for dust on the 
whole sky. Similarly, the covariance matrix is given by 
(c.f., Equation (fTTjl) 

C/ \ m „tcnsor , „scalar 

I j{r,s,a l ) = rc IJ + sc 



I J 

Ni,u + N 2 ,u 



(1 — ckd — "|)(1 - <*d - Q!g) ' ^ ^ 

where C/j, cjj and Nu denote a block of matrices for 
pixels within regions / and J. 

The free parameters in the maximization are r, s, ajj, 
and aig (i=l. . . 12(48)). In principle we wish to maximize 
the full likelihood function with respect to these param- 
eters; however, in practice, this process is too time con- 
suming to do brute-force, as varying each of these 15(51) 
parameters requires re-inverting a 4518 x 4518 matrix. 
Therefore, we make one approximation: we fix qd and 
a§ in the covariance matrix (Equation (|23p ) at the best- 
fit values, q;q and ctg°. This is a good approximation as 
long as the noise term is sub-dominant compared to the 
dominant scalar E'-modc signal, which is always the case 
for our low-noise configurationQ With this approxima- 
tion, 



£(r, s, oii) oc 



exp 



\x(cti) C(r,s;o$) x{at) 



(24) 



can be maximized with respect to r, s and on where i runs 
from 1 (dust) to 13(49) (s ynchrotron fo r 2 to 13(49)). Wc 
use the MINUIT package (|Jameslll988D for the maximiza- 
tion. 

In the fourth and fifth columns of Table [21 we show 
the recovered values of r for the synchrotron- only cases, 
in order to see if dividing the sky into 48 regions helps to 

6 As we have shown in Section 15.11 the foreground amplitudes 
and the dominant scalar E modes are covariant. Therefore, in 
order to find the best-fit as without running the full likelihood, we 
had to "cheat" and measure as in maps that do not contain the 
CMB signal or noise. Of course, we cannot do this in real life and 
thus we will have to come up with an efficient numerical algorithm 
for maximizing the full likelihood without this approximation. We 
believe that this is doable, so this will not be a limiting factor for 
our method. 
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Fig. 9. — Distribution of r and s obtained from 68 realizations 
input = 0.003 using Method I. 
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0.0000 
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r —r -4- 001 ft 
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Dust and Synchrotron (Method II) 
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input 



TABLE 3 

Recovered r from Synchrotron and Dust cleaning 



r- ^ a 
' input 


mean(r) b 


std(r) c 


mean(r) d 


std(r) c 


0.001 


0.0027 


0.0005 


0.0016 


0.0006 


0.003 


0.0050 


0.0008 


0.0038 


0.0009 


0.010 


0.0121 


0.0013 


0.0113 


0.0015 


0.030 


0.0326 


0.0021 






0.100 


0.1029 


0.0053 







a Input values of the scalar-to-tcnsor ratio for simulations 
(64(128) realizations for each r; nput for Method 1(H)). 
b Mean of the recovered maximum likelihood values of r 
for Method I. 

c Standard deviation of the recovered maximum likeli- 
hood values of r for Method I. 

d Mean of the recovered maximum likelihood values of r 
for Method II. 

c Standard deviation of the recovered maximum likeli- 
hood values of r for Method II. 

reduce the bias in r that we have just seen. We find 
that the bias has reduced, but not by much: Ar = 
0.00f8 -» 0.0014, 0.0019 0.0016, and 0.0020 -> 0.0015 
for r input = 0.001, 0.003, and 0.01, respectively. This is 
probably due to the division not being tailored to match 
the distribution of synchrotron emission. While we keep 
this simple division and do not pursue a more complex 
division in this paper, we shall come back to this issue in 
the future work. 

5.3. Recovering r 

Now, we recover r from the full dust-plus-synchrotron 
cases. In Figure [9l we show the distribution of r and s 
for all of 68 realizations that we have run with n npu t = 
0.003 using Method I. In Table[3j we show the recovered 
values of r in the second and fourth columns. Comparing 



them to the input values, 



input : 



in the first column, we 



conclude that our Method I recovers r with a foreground- 
induced bias of Ar ps 0.002, which is consistent with the 
bias we have just seen from the synchrotron-only cases. 



Fig. 10. — Recovered values of r (mean(r)) and error bars 
(std(r)) as a function of r; nput . The green points with error 
bars show the recovered r from the dust-only results (Table H): 
the red (cyan) points with error bars show the synchrotron- plus- 
dust results using Method I(II) (Table [3j A systematic bias of 
Ar 0.002(0.0006) is seen for the synchrotron-plus-dust results 
using Method I(II), which can be described by r recovere( j = ri npu t + 
(0.0018 ± 0.0004)((0.0006 ± 0.0004)) (red(cyan) dotted line). We 
do not detect an offset for the dust-only results: (0.0000 ± 0.0003) 
(green dashed line) 



With Method II we recover r with a much smaller bias 
of Ar w 0.0006. We visualize our results in Figure [TOT 

The bias in r is important, but the uncertainty in the 
recovered r is equally important. In Figure SI we have 
shown the predicted fractional errors on the determina- 
tion of r, ovAmput, for idealistic full-sky, foreground- 
free cases. How would they look when synchrotron and 
dust are included and cleaned with our method? In Fig- 
ure 1 1 1 1 we show the same figure but with the predic- 
tions made for the simulated CMB-plus-noise maps at 
^Vgide = 16 as described in Section [21 and scaled to the 
P06 mask. We also show er r /nnput where ay is extracted 
from the simulations. First, when both the foreground 
and lensing noise are ignored, the simulation and the an- 
alytical prediction are in a good agreement (to within 
10%) for r input = 0.001 (see the star symbol in Fig- 
ure [TT]). When the foreground is included, however, the 
error increases. For r; nput = 0.01, the foreground clean- 
ing increases the error by about 60%. We see larger dis- 
crepancies between the foreground-free predictions and 
the foreground-cleaned results for r input = 0.001: the 
foreground-cleaned error is a factor of two larger than 
the foreground-free prediction; thus, the increase in the 
error due to foreground cleaning can be substantial when 
r is as small as 10 -3 . 

Further optimizations could be done, given the details 
of a given experiment, and we intend to explore this issue 
within the context of some specific experimental designs. 
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Fig. 11. — The same as Figure[4] but the foreground-free predic- 
tions are made for simulated maps at N a i^ e = 16 (as described in 
Section [2]l and divided by ^/fsky where f a k y = 0.735 for the P06 
mask at N B ^ C = 16. The filled circles use the error bars, cr r , derived 
from our simulations with synchrotron and dust foreground clean- 
ing (denoted as "std(r)" in Table [31. The denominators, rj nput , arc 
not the values recovered from the simulations, but the input values 
fixed at either 0.01 or 0.001. Each filled circle (foreground-cleaned 
result from simulations) should be compared to the corresponding 
line with the same color (foreground- free, analytical prediction). 
The star shows the foreground- free and lensing-frec result from the 
simulation with r; nput = 0.001, which agrees with the analytical 
prediction to within 10%. 

Another improvement can be made by using N s id e = 32 
(or 64) map, so that the Nyquist frequency is close to (or 
beyond) the second bump of the B-mode spectrum and 
more information is used. Such analysis, however, would 
take 2 6 = 64 (4 6 = 4096) times more computation. 

6. CONCLUSIONS 

In this paper, we have studied the pixel-based fore- 
ground cleaning method within the context of a next- 
generation, low-noise CMB polarization satellite. This 
method was o riginally applied t o polarized data by th e 
WMAP team ([Page et alJl2007t IGold et al.ll2TM [2010h . 



and further investigated by lEfstathiou et al.l (|2009[) in 
the context of Planck. 

Despite the simplicity of the method (namely, we have 
maps at 3 different frequencies, two of which are used 
for removing the synchrotron and dust emission) , we are 
able to recover the input tensor-to-scalar ratio with only 
a small bias, Ar w 0.002(< 0.001) for the P06(extended) 
mask, which is dominated by the residual synchrotron 
emission. Further improvements should be straightfor- 
ward: one can tune the Galactic mask, and divide the 
synchrotron fitting regions according to the actual distri- 
bution of the synchrotron spectral index in the Galaxy 
(rather than using the regular division shown in Fig- 
ure [8]). One may also increase the number of frequencies 
for measuring the spatial distribution of the synchrotron 
spectral index, provided that we have enough space on 
the focal plane. These will be investigated in the context 
of specific experimental designs such as LiteBIRrQ, and 
presented elsewhere. 

Our study suggests that a detection of the primordial 
B-mode polarization at the level of r ~ 10~ 3 should be 
possible with carefully optimized mask and a regions. 
Note that our statistical error and systematic bias be- 
comes comparable with f s k y — 50% mask case. How- 
ever, let us mention one important caveat in our analysis. 
While our knowledge of the distribution and properties 
of the polarized synchrotron is fairly secure thanks to the 
WMAP data, our knowledge of the polarized dust emis- 
sion, especially the spatial variation of the dust spectral 
index, is still highly limited. Therefore, the estimated 
bias in r that we have presented in this paper cannot 
be too accurate. Fortunately, Planck will soon provide 
us with maps of the polarized dust emission with the 
unprecedented sensitivity; thus, we intend to revisit this 
issue once the Planck data become available. 
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APPENDIX 
A. SIGNAL COVARIANCE MATRIX 

Given power spectra, cf B and cf E , components of the signal covariance matrix for Q and U can be computed 
analytically. We have 



c(n n') = I C QQ(">') c QU (n,h') 
' \cuQ(n,n') cuu(n,n') J ' 



where 



c QQ (n, n') = £ cf V ]T W lm (n)Wr m (fi') + ]T cf % 2 £ X lm (h)Xf m (n') 

I m I m 

c QU (n, h') = J2 <f V Y,\-Wim{n)Xt m (ti)] + £ cf B wf £ X lm (h)w; m (n>) 

I m I m 

CUQ (h, h')=Y, c f Ew2 i £[-^ m (W ro (n')] + £ cf B wf ]T w lm (n)x; m (h') 

I m I m 

7 Light satellite for the studies of B-mode polarization I htt ■ / / mb kek ' — /lit bird l 
and Inflation from cosmic background Radiation Detection; 
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Cuu (n, A') = J2 c f V E Xim(n)XUn') + E c f V E W lm {n)Wt m {n') 

I m I m 

and 

W lm (n) = (-l)[ 2 y ;m (n) +_ 2 Y lm (n)}/2, 
X lm (n) = (-i)[ 2 Y lm (n) -_ 2 Y lm (n)]/2. 

We have assumed that E modes and B modes are uncorrelated. Here, wi is a smoothing function which includes an 
experimental beam, a pixel window function, and any other smoothing applied to maps. 

B. EXTENDED MASK 

The resolution 4 (r4) mask is extended from the P06 mask by setting the threshold foreground polarization intensity 
values at 60 and 240 GHz above which the pixels are masked. The intensity of the pixel i in the resolution 7 map is 
defined as 

Pii") = y/QlW + UH") (Bl) 

where Q and U are the sum of synchrotron and dust: 

[Qi, Ui](y) = [<2i, B ynch, Ui tSyndl }{v) + [Q. t ,dust, U iidnst ](i>) (B2) 

using PSM (See Eqs. G3) and CLl). 
An r4 pixel is masked if 

1. median of P^(240) pixels in the r4 pixel exceeds Threshold I, or 

2. maximum of P^(240) in the r4 pixel exceeds Threshold II, or 

3. median of P;(60) pixels in the r4 pixel exceeds Threshold III, or 

4. maximum of Pi (60) in the r4 pixel exceeds Threshold IV. 

Keeping f s k y = 50%, the values of the four thresholds are determined by minimizing the total foreground intensity in 
the residual map; 

P ros (mask)= \f®Ui + U L,i (B3) 

i^mask 

where 

[Qres,i, CU,i] = [QM(^)-aD[Qi,Ui](24D)-as[Qi,Ui](60). (B4) 
qd and as are given in the usual way by solving 

r=°.0' = AS) (B5) 

where 

X 2 = [QreB,<, ^ro S ,i] T [Qres,i, ^res.i], (Qres.i = ^res.t = 0, i G mask) (B6) 

The median and max. thresholds for the 240(60) GHz map determined this way are 19.2(1.42) and 38.4(2.11)/iA'. 
I.e., Threshold I = 19.2, II = 38.4, III = 1.42, and IV = 2.11 /iK. Note that we have defined an extended mask 
by using PSM maps without CMB or noise. In practice, both contributions would add noise spikes to the mask 
which need to be carefully examined. The noise contribution should be quite small given that we consider a low-noise 
(2 /iK arcmin) experiment in this paper. The CMB contribution can be removed by taking the difference between 
different channels and defining the threshold values on the difference maps (in the same way that the WMAP team has 
created temperature masks). However, given f s ky, Pes (mask) has a very broad bottom as a function of the thresholds. 
At the bottom, the shape of the mask is stable and our results are insensitive to the choice of the threshold values or 
the algorithm. 
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